% close all;
clear all;
load('colors.mat');

%% parameters
rootdir = ['..' filesep '191015_GeraldDoreenBob_Bats1'];
sampleNames = {};

pixSizeReal = 0.34; % um per pix

areaRatioFilterLevel = 0.7; % 0...1 how many data points inside the channel have to be valid in terms of area ratio
areaRatioMax = 1.2;

%% preamble
if rootdir(end) ~= filesep
    rootdir(end+1) = filesep;
end
if isempty(sampleNames)
    listDir = dir(rootdir);
    sampleNames = cell(length(listDir)-2,1);
    for k = 3:length(listDir)
        if exist([rootdir listDir(k).name filesep 'Online' filesep],'dir')
            sampleNames{k-2} = listDir(k).name;
        end
    end
    sampleNames = sampleNames(~cellfun('isempty',sampleNames));
end

%% choose samples
fprintf('Index\tName\n');
mask = true(1,length(sampleNames));
for k = 1:length(sampleNames)
    fprintf('%d\t%s\n',k,sampleNames{k});
end
idx = input('Choose samples: ');
if not(isempty(idx))
    sampleNames = sampleNames(idx);
end

%% extract cell traces
for idxSampleName = 1:length(sampleNames)
sampleName = sampleNames{idxSampleName};
filedir = [rootdir sampleName filesep 'Online' filesep];
[measNums,frameNums] = funcFindMeas(filedir);

pixSizeRec = pixSizeReal;%0.34;    % um per pix
exportMatName = '_data_export';%_contour';

firstImageBeforeChannelBegin = 1;
lastImageAfterChannelEnd     = 1;

createResultSubDir = 1;
writeStatImages = 1;
writeImages = 0;
enableFlag = 0;

cellNumberMat = zeros(length(measNums), 5);
cellNumberMat(:,1) = measNums(:);

diary([rootdir sampleName filesep 'analysis.log']);
fprintf('\n');
disp(datetime('now','TimeZone','Europe/Zurich','Format','d-MMM-y HH:mm:ss Z'));
fprintf('%s = %s\n','rootdir',rootdir);
fprintf('%s = %s\n','sampleName',sampleName);
fprintf('%s = ','measNums');fprintf('%d ',measNums);fprintf('\n');
fprintf('%s = ','frameNums');fprintf('%d ',frameNums);fprintf('\n');
fprintf('%s = %g\n','pixSizeReal',pixSizeReal);
fprintf('%s = %g\n','pixSizeRec',pixSizeRec);
fprintf('%s = %s\n','exportMatName',exportMatName);
fprintf('%s = %d\n','firstImageBeforeChannelBegin',firstImageBeforeChannelBegin);
fprintf('%s = %d\n','lastImageAfterChannelEnd',lastImageAfterChannelEnd);
fprintf('%s = %d\n','createResultSubDir',createResultSubDir);
fprintf('%s = %d\n','writeStatImages',writeStatImages);
fprintf('%s = %d\n','writeImages',writeImages);
fprintf('%s = %d\n','enableFlag',enableFlag);
fprintf('%s = %g\n','areaRatioFilterLevel',areaRatioFilterLevel);
fprintf('%s = %g\n','areaRatioMax',areaRatioMax);

for k = 1:5
    figure(k);clf;
end

for measNumN = 1:length(measNums)
measNum = measNums(measNumN);
fps = frameNums(measNumN);

%% initialization
file = dir([filedir sprintf('M%i_*.tdms', measNum)]);
filename = file.name;
dt = 1/fps;

filesep_idx = strfind(rootdir, filesep);
underscore_idx = strfind(rootdir(filesep_idx(end-1)+1:end), '_');
if underscore_idx(1) <= 9
    dateLength = underscore_idx(1)-1;
else
    tmp = rootdir(filesep_idx(end-1)+1:end);
    if (0 <= (tmp(7)-48)) && ((tmp(7)-48) <= 9) && (0 <= (tmp(8)-48)) && ((tmp(8)-48) <= 9)
        dateLength = 8;
    else
        dateLength = 6;
    end
end
date = rootdir(filesep_idx(end-1)+1:filesep_idx(end-1)+1+dateLength-1);
measName = rootdir(filesep_idx(end-1)+1+dateLength+1:end);
clear filesep_idx

fileChPos = dir([filedir sprintf('M%i_channelPos.mat',measNum)]);
if isempty(fileChPos)
    fileChPos = dir([filedir 'channelPos.mat']); % old format
    load([filedir fileChPos.name]);
else
    load([filedir fileChPos.name]);
    if exist('chPos','var')
        channelPos = chPos;
    else
%         channelPos already exists
    end
end

%% import data from tdms file
tdmsdata = convertTDMS(0,[filedir filename]);
datastruct = tdmsdata.Data.MeasuredData;
vFrame = datastruct(3).Data';
vXPos = datastruct(4).Data' *pixSizeReal/pixSizeRec;
vYPos = datastruct(5).Data' *pixSizeReal/pixSizeRec;
vAx1 = datastruct(6).Data' *pixSizeReal/pixSizeRec;
vAx2 = datastruct(7).Data' *pixSizeReal/pixSizeRec;
vCircularity = datastruct(8).Data';
vArea = datastruct(9).Data';
vRawArea = datastruct(10).Data';
vXPosAcc = datastruct(11).Data' *pixSizeReal/pixSizeRec;
vCellNum = datastruct(12).Data';

img2tdmsOffset = -1;

datalength = length(vFrame);
clear tdmsdata datastruct;

%% find cell sequences
vCellIdx = double(diff(vCellNum) > 0) .* (1:length(vCellNum)-1);
vCellIdx(vCellIdx == 0) = [];
vCellLength = diff([0 vCellIdx length(vCellNum)]);  % length of cell sequences
vCellIdx = [1 vCellIdx+1];  % start index for cell sequence

set(0,'CurrentFigure',1);
plot(1:length(vXPosAcc),vXPosAcc, 1:length(vXPosAcc),repmat(channelPos(1)*pixSizeReal,1,length(vXPosAcc)),'r', 1:length(vXPosAcc),repmat(channelPos(2)*pixSizeReal,1,length(vXPosAcc)),'r');
xlabel('no. of datapoint');ylabel('accumulated x position in um');

set(0,'CurrentFigure',2);
hist(vCellLength,20);
ylabel('count');xlabel('lenght of cell sequence');

% CellLenghtTh = input('Threshold = ');
CellLenghtTh = 10;
vCellIdx(vCellLength < CellLenghtTh) = [];
vCellLength(vCellLength < CellLenghtTh) = [];

cellNumberMat(measNumN,2) = length(vCellIdx);
fprintf('min. %i cells per sequence: %i cell sequences\n', CellLenghtTh, length(vCellIdx));

%% filter cell sequences for conditions regarding first and last cell position
if( firstImageBeforeChannelBegin || lastImageAfterChannelEnd )
    filterMask = false(1, length(vCellIdx));
    for cellNumPtr = 1:length(vCellIdx)
        mask = (vCellIdx(cellNumPtr):vCellIdx(cellNumPtr)+vCellLength(cellNumPtr)-1);
        vXPosAcc_msk = vXPosAcc(mask);

        if( firstImageBeforeChannelBegin && (vXPosAcc_msk(1) > channelPos(1)*pixSizeReal) )
            filterMask(cellNumPtr) = true;
        end
        if( lastImageAfterChannelEnd && (vXPosAcc_msk(end) < channelPos(2)*pixSizeReal) )
            filterMask(cellNumPtr) = true;
        end
    end
    vCellIdx(filterMask) = [];
    vCellLength(filterMask) = [];
    cellNumberMat(measNumN,3) = length(vCellIdx);
    fprintf('cell traces covers the whole channel: %i cell sequences\n', length(vCellIdx));
end

%% calculate time
vTime = (vFrame - vFrame(1)) * dt * 1e3;    % in ms

%% load video
% if writeImages
    fprintf('load video\n');
    fileAvi = dir([filedir sprintf('M%i_*imaq.avi', measNum)]);
    vr = VideoReader([filedir fileAvi.name]);
    numFrames = get(vr,'NumberOfFrames');
    frame = read(vr,1);
    clear fileAvi;
% end

%% load contours
fprintf('load contours\n');
[contourData, contourDataIdx] = readCellContours( filedir, measNum, numFrames+10 );

fprintf('... and calculate Lambda ratio, volume, surface area\n');
vFramePos = ones(2,length(vFrame));
vLambdaRatio = zeros(1,length(vFrame));
vVolume = zeros(1,length(vFrame));
vSurfaceArea = zeros(1,length(vFrame));
for l = 1:length(vFrame)
    framePos = contourDataIdx(contourDataIdx(:,1) == vFrame(l),2:3);
    if ~isempty(framePos)
        framePos = framePos(1,:);
        vFramePos(:,l) = framePos(:);
        [ T, area, Ix, Iy, xb, yb ] = calcMomentOfAreaTensor(contourData(framePos(1):framePos(1)+framePos(2)-1,:));
        [v,d] = eig(T);
        eigenvals = diag(d);
        mask = abs(v(1,:)) > abs(v(2,:));
        if sum(mask) == 0
            mask = [true false];
        end
        e = [eigenvals(mask), eigenvals(not(mask))];
        vLambdaRatio(l) = e(1) / e(2);
        [vVolume(l), vSurfaceArea(l)] = calcVolumeSurfacearea(contourData(framePos(1):framePos(1)+framePos(2)-1,:));
    end
end

%% results
if createResultSubDir
    filedirImg = [filedir sprintf('M%i_eval_result', measNum) filesep];
    if ~exist(filedirImg,'dir')
        mkdir(filedirImg);
    end
else
    filedirImg = filedir;
end

if writeStatImages
    print(1,[filedirImg sprintf('M%i_eval_vXPosAcc', measNum)], '-dpng');
    print(2,[filedirImg sprintf('M%i_eval_histCellSequences', measNum)], '-dpng');
end

set(0,'CurrentFigure',3);
figPos = get(gcf, 'Position');
set(gcf, 'Position', [figPos(1) figPos(2) 1120 630]);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperOrientation', 'portrait');
set(gcf, 'PaperPositionMode', 'manual');
res = 100;
set(gcf, 'PaperSize', [1120/res*2.54 630/res*2.54]);
set(gcf, 'PaperPosition', [0 0 1120/res*2.54 630/res*2.54]);

%% filter for area ratio
maskExcludeCells = false(1,length(vCellIdx));
for cellNumPtr = 1:length(vCellIdx)
    mask = (vCellIdx(cellNumPtr):vCellIdx(cellNumPtr)+vCellLength(cellNumPtr)-1);
    vAreaRatioMask = vArea(mask)./vRawArea(mask);
    maskInsideChannel = (channelPos(1)*pixSizeReal <= vXPosAcc(mask)) & (vXPosAcc(mask) <= channelPos(2)*pixSizeReal);
    vAreaRatioMask = vAreaRatioMask(maskInsideChannel);
    % less than x% valid points
    if sum(double(vAreaRatioMask <= areaRatioMax))/length(vAreaRatioMask) < areaRatioFilterLevel
        maskExcludeCells(cellNumPtr) = true;
    end
end
fprintf('filter for area ratio: excluded cells: %u\n', sum(double(maskExcludeCells)));
vCellIdx(maskExcludeCells) = [];
vCellLength(maskExcludeCells) = [];
cellNumberMat(measNumN,4) = length(vCellIdx);

%% 
cellData = repmat(struct('number',[], 'videoPointers',[], 'data',[]),1,length(vCellIdx));
fprintf('write results (%d cells):\n', length(vCellIdx));
divCellNumPtr = length(vCellIdx)/10;
fprintf('----------\n');
chCoverageEnd = zeros(1,length(vCellIdx));
velocities = zeros(1,length(vCellIdx));
velTimes = zeros(1,length(vCellIdx));
vFlagExclude = false(1,length(vCellIdx));
for cellNumPtr = 1:length(vCellIdx)
    mask = (vCellIdx(cellNumPtr):vCellIdx(cellNumPtr)+vCellLength(cellNumPtr)-1);
    cellNum = vCellNum(vCellIdx(cellNumPtr)+round(0.5*(vCellLength(cellNumPtr)-1)));
    
    videoPointers = [max(vCellIdx(cellNumPtr)+img2tdmsOffset,1) min(vCellIdx(cellNumPtr)+vCellLength(cellNumPtr)-1+img2tdmsOffset,numFrames)];
    
    vXPosAcc_msk = vXPosAcc(mask);
    chCoverageEnd(cellNumPtr) = (vXPosAcc_msk(end)/pixSizeReal - channelPos(1)) / (channelPos(2) - channelPos(1));
    maskInsideChannel = (channelPos(1)*pixSizeReal <= vXPosAcc_msk) & (vXPosAcc_msk <= channelPos(2)*pixSizeReal);
    vTime_msk = vTime(mask);
    medianVelocity = nanmedian(gradient(vXPosAcc_msk(maskInsideChannel), vTime_msk(maskInsideChannel))*0.1);
    velocities(cellNumPtr) = medianVelocity;
    velTimes(cellNumPtr) = nanmedian(vTime_msk);
    
    vFramePos_msk = vFramePos(:,mask);
    contourData_msk = contourData(vFramePos_msk(1,1):vFramePos_msk(1,end)+vFramePos_msk(2,end)-1,:);
    contourDataIdx_msk = [(vFramePos_msk(1,:)-vFramePos_msk(1,1)+1)', (vFramePos_msk(2,:))'];
    
    cellDataData = struct('time',vTime(mask), 'frame',vFrame(mask), 'xPos',vXPos(mask), 'yPos',vYPos(mask), 'ax1',vAx1(mask), 'ax2',vAx2(mask), 'circularity',vCircularity(mask), 'area',vArea(mask)*pixSizeReal^2, 'areaRaw',vRawArea(mask)*pixSizeReal^2, 'xPosAcc',vXPosAcc_msk, 'lambdaRatio',vLambdaRatio(mask), 'areaRatio',vArea(mask)./vRawArea(mask), 'volume',vVolume(mask)*pixSizeReal^3, 'surfaceArea',vSurfaceArea(mask)*pixSizeReal^2,...
        'medianArea',median(vArea(mask)*pixSizeReal^2), 'medianAreaRaw',median(vRawArea(mask)*pixSizeReal^2), 'medianVelocity',medianVelocity, ...
        'contourIndex',contourDataIdx_msk, 'contourData',contourData_msk);
    cellData(cellNumPtr).number = cellNum;
    cellData(cellNumPtr).videoPointers = videoPointers;
    cellData(cellNumPtr).data = cellDataData;
    
    if writeImages
        % 6 images
        xCoord = vXPosAcc_msk - channelPos(1)*pixSizeReal;
        endCoord = diff(channelPos)*pixSizeReal;
        % first frame
        idxFrame = 1;   frames{1} = read(vr,videoPointers(1)+idxFrame-1);cntrs{1} = contourData_msk(contourDataIdx_msk(idxFrame,1):contourDataIdx_msk(idxFrame,1)+contourDataIdx_msk(idxFrame,2)-1,:);
        % frame @ inlet
        [~,idxFrame] = min(abs(xCoord - 0*endCoord));
        frames{2} = read(vr,videoPointers(1)+idxFrame-1); cntrs{2} = contourData_msk(contourDataIdx_msk(idxFrame,1):contourDataIdx_msk(idxFrame,1)+contourDataIdx_msk(idxFrame,2)-1,:);
        % @ 50% in channel
        [~,idxFrame] = min(abs(xCoord - 0.5*endCoord));     frames{3} = read(vr,videoPointers(1)+idxFrame-1);cntrs{3} = contourData_msk(contourDataIdx_msk(idxFrame,1):contourDataIdx_msk(idxFrame,1)+contourDataIdx_msk(idxFrame,2)-1,:);
        % @ 95% in channel
        [~,idxFrame] = min(abs(xCoord - 0.95*endCoord));    frames{4} = read(vr,videoPointers(1)+idxFrame-1);cntrs{4} = contourData_msk(contourDataIdx_msk(idxFrame,1):contourDataIdx_msk(idxFrame,1)+contourDataIdx_msk(idxFrame,2)-1,:);
        % last frame before passing outlet
        [~,idxFrame] = max(xCoord .* (xCoord <= endCoord));     frames{5} = read(vr,videoPointers(1)+idxFrame-1);cntrs{5} = contourData_msk(contourDataIdx_msk(idxFrame,1):contourDataIdx_msk(idxFrame,1)+contourDataIdx_msk(idxFrame,2)-1,:);
        % last frame
        frames{6} = read(vr,videoPointers(end));cntrs{6} = contourData_msk(contourDataIdx_msk(end,1):contourDataIdx_msk(end,1)+contourDataIdx_msk(end,2)-1,:);

        showResultFigure6( 3, cellData(cellNumPtr), measNum, sampleName, frames, cntrs, pixSizeReal, channelPos, areaRatioMax, [color1;color2;color3;color4] )
        print( 3,'-dpng','-painters','-r120',[filedirImg sprintf('M%i_cell%04i', measNum, cellNum)]);
        clear frames;
    end
    
    if enableFlag
        if input('take it? y/n: ','s') ~= 'y'
            vFlagExclude(cellNumPtr) = true;
        end
    else
        fprintf(repmat('-',1,round(cellNumPtr/divCellNumPtr) - round((cellNumPtr-1)/divCellNumPtr)));
    end
end
while sum(double(vFlagExclude)) > 0
    idx = find(vFlagExclude,1,'first');
    cellData(idx) = [];
    vFlagExclude(idx) = [];
end
fprintf('\n');
measurement = struct('number',measNum, 'fps',fps, 'pixelSize',pixSizeReal, 'cellData',cellData);
channelPositions = struct('pix',channelPos, 'um',channelPos*pixSizeReal);
sample = struct('measName',measName, 'name',sampleName, 'date',date, 'channelPositions',channelPositions, 'measData',measurement);

set(0,'CurrentFigure',4);
hist(chCoverageEnd);xlabel('end point in % of channel length');xlim([0.7 1.5]);
if writeStatImages
    print(4,[filedirImg sprintf('M%i_eval_chCoverageEnd', measNum)], '-dpng');
end

set(0,'CurrentFigure',5);
plot(velTimes,velocities,'.-');xlabel('time in ms');ylabel('median velocity in cm/s');
if writeStatImages
    print(5,[filedirImg sprintf('M%i_eval_velocity', measNum)], '-dpng');
end

%% export result file

save([filedir sprintf('M%i', measNum) exportMatName],'sample');
fprintf('%d cells written\n', length(sample.measData.cellData));
cellNumberMat(measNumN,5) = length(sample.measData.cellData);

end
%%
fprintf('measurement: %s\n',filedir);
fprintf('cell counts:\n');
fprintf('M#\tTh\twholeChannel\tareaRatio\twritten\n');
fprintf('%i\t%i\t%i\t%i\t%i\n', cellNumberMat');

diary off;

end

return
%% functions
function [measNums,frameNums] = funcFindMeas(fileDir)
fileList = dir([fileDir 'M*_camera.ini']);
measNums  = NaN(1,length(fileList));
frameNums = NaN(1,length(fileList));
strKeyFrame = 'Frame Rate = ';
for idxFile = 1:length(fileList)
    fid = fopen([fileDir fileList(idxFile).name],'r');
    frameNums(idxFile) = getFrame(fid, strKeyFrame);
    idx1 = strfind(fileList(idxFile).name, 'M');
    idx2 = strfind(fileList(idxFile).name, '_');
    measNums(idxFile) = str2double(fileList(idxFile).name(idx1(1)+1:idx2(1)-1));
    fclose(fid);
end
end

function frame = getFrame(fid, strKeyFrame)
while ~feof(fid)
    strLine = fgets(fid);
    index = strfind(strLine, strKeyFrame);
    if ~isempty(index)
        frame = str2double(strLine(index + length(strKeyFrame):end));
        return;
    end
end
end

